# TODO: Add comment
# 
# Author: yaping
###############################################################################

library(gridExtra)
library(ggplot2)

makeScatterHistplot<-function(x, name="test.scatterHisplot.png", loc=NULL, region=c(7,9),xlimit=NULL,ylimit=NULL,xlabel=NULL,ylabel=NULL, ...){
	#env <- new.env(parent = baseenv())
	png(name, height=600, width=600)
	#layout(rbind(c(2,1),c(0,3)), widths =c(1,2), heights =c(2,1), respect = FALSE)
	overlapped<-NULL
	overlapped_name<-NULL
	if(!is.null(loc)){		
		overlapped_name<-ifelse(countOverlaps(GRanges(seqnames=x$chr,ranges=IRanges(start=x$start,end=x$end),strand="*"),loc)>0,as.character(x$name),"")
		overlapped_name<-overlapped_name[!(overlapped_name %in% "")]
		overlapped<-data.frame(a=x[overlapped_name,][,region[1]],b=x[overlapped_name,][,region[2]])
	}else{
		overlapped<-data.frame(a=x[,1],b=x[,2])
	}	
	scatter <- ggplot(overlapped,aes(a, b))+geom_point(alpha=0.3)+ xlab(xlabel)+ylab(ylabel)+ xlim(xlimit[1], xlimit[2])+ ylim(ylimit[1], ylimit[2])
	hist_right <- ggplot(overlapped,aes(b))+geom_histogram()+coord_flip()+scale_x_continuous(breaks = ylimit[1]:ylimit[2],limits=c(ylimit[1], ylimit[2]),expand = c(.05,.05))+xlab(ylabel)
	hist_top <- ggplot(overlapped,aes(a))+geom_histogram()+scale_x_continuous(breaks = xlimit[1]:xlimit[2],limits=c(xlimit[1], xlimit[2]),expand = c(.05,.05))+ xlab(xlabel)
	empty <- ggplot()+geom_point(aes(1,1), colour="white")+
			theme(legend.position = "none", axis.title.x = element_blank(),axis.title.y = element_blank(),axis.text.x = element_blank(),axis.text.y = element_blank(),axis.line=element_blank(),panel.background = element_blank(),panel.grid.major = element_blank(),
					panel.grid.minor = element_blank(),panel.border = element_blank())
	
	#grid.arrange(hist_top, empty, scatter, hist_right, ncol=2, nrow=2, widths=c(2, 1), heights=c(1, 2))
	do.call(grid.arrange, list(hist_top, empty, scatter, hist_right, ncol=2, nrow=2, widths=c(2, 1), heights=c(1, 2)))
	dev.off()
}

makeScatterHistplot<-function(x, matrix2=NULL,name="test.scatterHisplot.png", loc=NULL, region=c(7,9),xlimit=NULL,ylimit=NULL,xlabel=NULL,ylabel=NULL, ...){
	#env <- new.env(parent = baseenv())
	png(name, height=600, width=600)
	#layout(rbind(c(2,1),c(0,3)), widths =c(1,2), heights =c(2,1), respect = FALSE)
	overlapped<-NULL
	overlapped_name<-NULL
	overlapped_2<-NULL
	if(!is.null(loc)){		
		overlapped_name<-ifelse(countOverlaps(GRanges(seqnames=x$chr,ranges=IRanges(start=x$start,end=x$end),strand="*"),loc)>0,as.character(x$name),"")
		overlapped_name<-overlapped_name[!(overlapped_name %in% "")]
		overlapped<-data.frame(a=x[overlapped_name,][,region[1]],b=x[overlapped_name,][,region[2]])
		if(!is.null(matrix2)){
			overlapped_2<-data.frame(c=matrix2[intersect(overlapped_name,rownames(matrix2)),][,region[1]],d=matrix2[intersect(overlapped_name,rownames(matrix2)),][,region[2]])
		}
	}else{
		overlapped<-data.frame(a=x[,1],b=x[,2])
		if(!is.null(matrix2)){
			overlapped_2<-data.frame(c=matrix2[,1],d=matrix2[,2])
		}
	}	
	scatter <- ggplot(overlapped,aes(a, b))+geom_point(alpha=0.3)+ xlab(xlabel)+ylab(ylabel)+ xlim(xlimit[1], xlimit[2])+ ylim(ylimit[1], ylimit[2])
	
	hist_right <- ggplot(overlapped,aes(b))+geom_histogram()+coord_flip()+scale_x_continuous(breaks = ylimit[1]:ylimit[2],limits=c(ylimit[1], ylimit[2]),expand = c(.05,.05))+xlab(ylabel)
	hist_top <- ggplot(overlapped,aes(a))+geom_histogram()+scale_x_continuous(breaks = xlimit[1]:xlimit[2],limits=c(xlimit[1], xlimit[2]),expand = c(.05,.05))+ xlab(xlabel)
	if(!is.null(overlapped_2)){
		scatter<-scatter+geom_point(alpha=0.8,data=overlapped_2,aes(c,d,colour = "red"))+theme(legend.position = "none")
		hist_right <- ggplot()+geom_histogram(data=overlapped,aes(b),fill = "black", alpha = 0.4)+geom_histogram(data=overlapped_2,aes(d),fill = "red", alpha = 0.4)+coord_flip()+scale_x_continuous(breaks = ylimit[1]:ylimit[2],limits=c(ylimit[1], ylimit[2]),expand = c(.05,.05))+xlab(ylabel)
		hist_top <- ggplot()+geom_histogram(data=overlapped,aes(a),fill = "black", alpha = 0.4)+geom_histogram(data=overlapped_2,aes(c),fill = "red", alpha = 0.4)+scale_x_continuous(breaks = xlimit[1]:xlimit[2],limits=c(xlimit[1], xlimit[2]),expand = c(.05,.05))+ xlab(xlabel)
	}
	empty <- ggplot()+geom_point(aes(1,1), colour="white")+
			theme(legend.position = "none", axis.title.x = element_blank(),axis.title.y = element_blank(),axis.text.x = element_blank(),axis.text.y = element_blank(),axis.line=element_blank(),panel.background = element_blank(),panel.grid.major = element_blank(),
					panel.grid.minor = element_blank(),panel.border = element_blank())
	
	#grid.arrange(hist_top, empty, scatter, hist_right, ncol=2, nrow=2, widths=c(2, 1), heights=c(1, 2))
	do.call(grid.arrange, list(hist_top, empty, scatter, hist_right, ncol=2, nrow=2, widths=c(2, 1), heights=c(1, 2)))
	dev.off()
}


#e.g. makeScatterHistplot(cbind(log(rpkm.1stExon[,1]),log(rpkm.1stExon[,2])),name="rep1_vs_rep2(log(RPKM) in 1st exon, HCT116).scatterplot.png",xlabel="rep1",ylabel="rep2",xlimit=c(-10,10),ylimit=c(-10,10))


makeMergedHistplot<-function(x, matrix2=NULL,name="test.mergedHisplot.png", loc=NULL, loc2=NULL, region=c(7,9),xlimit=NULL,xlabel=NULL){
	#env <- new.env(parent = baseenv())
	png(name, height=600, width=600)
	#layout(rbind(c(2,1),c(0,3)), widths =c(1,2), heights =c(2,1), respect = FALSE)
	overlapped<-NULL
	overlapped_name<-NULL
	overlapped_2<-NULL
	overlapped_loc2<-NULL
	if(!is.null(loc)){		
		overlapped_name<-ifelse(countOverlaps(GRanges(seqnames=x$chr,ranges=IRanges(start=x$start,end=x$end),strand="*"),loc)>0,as.character(x$name),"")
		overlapped_name<-overlapped_name[!(overlapped_name %in% "")]
		overlapped<-data.frame(a=x[overlapped_name,][,region[1]],b=x[overlapped_name,][,region[2]])
		if(!is.null(loc2)){
			
			#overlapped$b=x[overlapped_name_2,][,region[2]]
		}
		if(!is.null(matrix2)){
			overlapped_2<-data.frame(c=matrix2[intersect(overlapped_name,rownames(matrix2)),][,region[1]],d=matrix2[intersect(overlapped_name,rownames(matrix2)),][,region[2]])
			
		}else{
			if(!is.null(loc2)){
				overlapped_name_2<-ifelse(countOverlaps(GRanges(seqnames=x$chr,ranges=IRanges(start=x$start,end=x$end),strand="*"),loc2)>0,as.character(x$name),"")
				overlapped_name_2<-overlapped_name_2[!(overlapped_name_2 %in% "")]
				overlapped_loc2<-data.frame(a=x[overlapped_name_2,][,region[1]],b=x[overlapped_name_2,][,region[2]])
			}
		}
	}else{
		overlapped<-data.frame(a=x[,1],b=x[,2])
		if(!is.null(matrix2)){
			overlapped_2<-data.frame(c=matrix2[,1],d=matrix2[,2])
		}
	}	
	
	hist <- ggplot()+geom_histogram(data=overlapped,aes(a),fill = "black", alpha = 0.4)+geom_histogram(data=overlapped,aes(b),fill = "blue", alpha = 0.4)+scale_x_continuous(breaks = xlimit[1]:xlimit[2],limits=c(xlimit[1], xlimit[2]),expand = c(.05,.05))+ xlab(xlabel)
	if(!is.null(overlapped_2)){
		hist <- ggplot()+geom_histogram(data=overlapped,aes(a),fill = "black", alpha = 0.4)+geom_histogram(data=overlapped,aes(b),fill = "blue", alpha = 0.4) +geom_histogram(data=overlapped_2,aes(c),fill = "red", alpha = 0.4)+geom_histogram(data=overlapped_2,aes(d),fill = "purple", alpha = 0.4)+scale_x_continuous(breaks = xlimit[1]:xlimit[2],limits=c(xlimit[1], xlimit[2]),expand = c(.05,.05))+ xlab(xlabel)
	}else{
		if(!is.null(loc2)){
			hist <- ggplot()+geom_histogram(data=overlapped,aes(a),fill = "black", alpha = 0.4)+geom_histogram(data=overlapped_loc2,aes(a),fill = "blue", alpha = 0.4)+scale_x_continuous(breaks = xlimit[1]:xlimit[2],limits=c(xlimit[1], xlimit[2]),expand = c(.05,.05))+ xlab(xlabel)
		}
	}
	
	do.call(grid.arrange, list(hist, ncol=1, nrow=1, widths=1, heights=1))
	dev.off()
}

makeMergedDensityplot<-function(x, matrix2=NULL,name="test.mergedDensityplot.png", loc=NULL, loc2=NULL, region=c(7,9),xlimit=NULL,xlabel=NULL,ylimit=NULL){
	#env <- new.env(parent = baseenv())
	png(name, height=600, width=600)
	#layout(rbind(c(2,1),c(0,3)), widths =c(1,2), heights =c(2,1), respect = FALSE)
	overlapped<-NULL
	overlapped_name<-NULL
	overlapped_2<-NULL
	overlapped_loc2<-NULL
	if(!is.null(loc)){		
		overlapped_name<-ifelse(countOverlaps(GRanges(seqnames=x$chr,ranges=IRanges(start=x$start,end=x$end),strand="*"),loc)>0,as.character(x$name),"")
		overlapped_name<-overlapped_name[!(overlapped_name %in% "")]
		overlapped<-data.frame(a=x[overlapped_name,][,region[1]],b=x[overlapped_name,][,region[2]])
		if(!is.null(loc2)){
			
			#overlapped$b=x[overlapped_name_2,][,region[2]]
		}
		if(!is.null(matrix2)){
			overlapped_2<-data.frame(c=matrix2[intersect(overlapped_name,rownames(matrix2)),][,region[1]],d=matrix2[intersect(overlapped_name,rownames(matrix2)),][,region[2]])
			
		}else{
			if(!is.null(loc2)){
				overlapped_name_2<-ifelse(countOverlaps(GRanges(seqnames=x$chr,ranges=IRanges(start=x$start,end=x$end),strand="*"),loc2)>0,as.character(x$name),"")
				overlapped_name_2<-overlapped_name_2[!(overlapped_name_2 %in% "")]
				overlapped_loc2<-data.frame(a=x[overlapped_name_2,][,region[1]],b=x[overlapped_name_2,][,region[2]])
			}
		}
	}else{
		overlapped<-data.frame(a=x[,1],b=x[,2])
		if(!is.null(matrix2)){
			overlapped_2<-data.frame(c=matrix2[,1],d=matrix2[,2])
		}
	}	
	
	hist <- ggplot()+geom_density(data=overlapped,aes(a),color = "black")+geom_density(data=overlapped,aes(b),color = "blue")+scale_x_continuous(breaks = xlimit[1]:xlimit[2],limits=c(xlimit[1], xlimit[2]),expand = c(.05,.05))+ xlab(xlabel)+scale_y_continuous(breaks = seq(ylimit[1],ylimit[2],by=(ylimit[2]-ylimit[1])/5),limits=c(ylimit[1], ylimit[2]))
	if(!is.null(overlapped_2)){
		hist <- ggplot()+geom_density(data=overlapped,aes(a),color = "black")+geom_density(data=overlapped,aes(b),color = "blue") +geom_density(data=overlapped_2,aes(c),color = "red")+geom_density(data=overlapped_2,aes(d),color = "purple")+scale_x_continuous(breaks = xlimit[1]:xlimit[2],limits=c(xlimit[1], xlimit[2]),expand = c(.05,.05))+ xlab(xlabel)+scale_y_continuous(breaks = seq(ylimit[1],ylimit[2],by=(ylimit[2]-ylimit[1])/5),limits=c(ylimit[1], ylimit[2]))
	}else{
		if(!is.null(loc2)){
			hist <- ggplot()+geom_density(data=overlapped,aes(a),color = "black")+geom_density(data=overlapped_loc2,aes(a),color = "blue")+scale_x_continuous(breaks = xlimit[1]:xlimit[2],limits=c(xlimit[1], xlimit[2]),expand = c(.05,.05))+ xlab(xlabel)+scale_y_continuous(breaks = seq(ylimit[1],ylimit[2],by=(ylimit[2]-ylimit[1])/5),limits=c(ylimit[1], ylimit[2]))
		}
	}
	
	do.call(grid.arrange, list(hist, ncol=1, nrow=1, widths=1, heights=1))
	dev.off()
}


library(easyRNASeq)
library(edgeR)
library(DESeq)
library(NOISeq)
##read raw read count
count.table.allUniqExon <- easyRNASeq(filesDirectory="/home/yapingli/project/HCT116_DKO1/RNAseq/bam",filenames=c("HCT116_Fides_D28MLACXX_1_KEL656A236.male.hg19.nodups.uniq.bam","HCT116_Heather_D28MLACXX_1_WIT1251A69.male.hg19.nodups.uniq.bam","DKO1_Fides_D28MLACXX_1_KEL656A237.male.hg19.nodups.uniq.bam","DKO1_Heather_D28MLACXX_1_WIT1251A70.male.hg19.nodups.uniq.bam"),conditions=HCT116_conditions,readLength=50L,organism="Hsapiens",annotationFile="/home/yapingli/uec-gs1/data/genome_data/common_genomic_regions/hg19/knownGene.genePredToGtf.ExonOnly.uniq.transcriptId.hg19.gtf",count="exons",annotationMethod="gtf",gapped=TRUE)

##edgeR GLM test:
##e.g.resultByEdgeRGlmTest<-sigExonByEdgeRGlmTest(count.table.1stExon,anno=count.table.1stExon.GC.geneLen)
##e.g.resultByEdgeRGlmTestAllExon<-sigExonByEdgeRGlmTest(count.table.allUniqExon,anno=count.table.allUniqExon.GC.geneLen)
sigExonByEdgeRGlmTest<-function(count.table,group=factor(c(1,1,2,2)),filterNum=1,fdr=0.05, coef=2, anno=NULL){
	filter <- apply(count.table,1,function(x) sum(x)>=filterNum)
	common <- intersect(rownames(anno),rownames(count.table[filter,]))
	count.table.filter<-count.table[common,]
	anno.filter<-anno[common,]
	dge <- DGEList(counts=count.table.filter, group=group)
	design <- model.matrix(~group)
	dge <- estimateGLMCommonDisp(dge, design)
	dge <- estimateGLMTrendedDisp(dge, design)
	dge <- estimateGLMTagwiseDisp(dge, design)
	fit <- glmFit(dge, design)
	lrt <- glmLRT(fit, coef=coef)
	lrt.table<-lrt$table
	lrt.table<-cbind(lrt.table,p.adjust(lrt.table$PValue,method="BH"))
	up_regul.table<-lrt.table[lrt.table[,5]<fdr & lrt.table[,1]> 1,]
	down_regul.table<-lrt.table[lrt.table[,5]<fdr & lrt.table[,1] < (-1),]
	up_regul.genes<-anno.filter[rownames(up_regul.table),]
	up_regul.genes<-cbind(up_regul.genes[,1:3],up_regul.genes[,6],up_regul.genes[,5],up_regul.genes[,4])
	down_regul.genes<-anno.filter[rownames(down_regul.table),]
	down_regul.genes<-cbind(down_regul.genes[,1:3],down_regul.genes[,6],down_regul.genes[,5],down_regul.genes[,4])
	print(paste("Up regulated exon number is: ",dim(up_regul.table)[1]))
	print(paste("Down regulated exon number is: ",dim(down_regul.table)[1]))
	returnResult<-list(upRegulGenes=up_regul.genes, upRegulTables=up_regul.table, downRegulGenes=down_regul.genes, downRegulTables=down_regul.table)
	return(returnResult)
}

##DE-seq exact test:
##e.g.resultByDESeqExactTest<-sigExonByDESeqExactTest(count.table.1stExon,anno=count.table.1stExon.GC.geneLen)
##resultByDESeqExactTestAllExon<-sigExonByDESeqExactTest(count.table.allUniqExon,anno=count.table.allUniqExon.GC.geneLen)
sigExonByDESeqExactTest<-function(count.table,group=factor(c(1,1,2,2)),filterNum=1,fdr=0.05, anno=NULL){
	filter <- apply(count.table,1,function(x) sum(x)>=filterNum)
	common <- intersect(rownames(anno),rownames(count.table[filter,]))
	count.table.filter<-count.table[common,]
	anno.filter<-anno[common,]
	cds = newCountDataSet( count.table.filter, group )
	cds = estimateSizeFactors( cds )
	cds = estimateDispersions( cds )
	res = nbinomTest( cds, "1", "2" )
	up_regul.table<-res[res$foldChange > 2 & res$padj < fdr,]
	down_regul.table<-res[res$foldChange < 0.5 & res$padj < fdr,]
	up_regul.genes<-anno.filter[up_regul.table$id,]
	up_regul.genes<-cbind(up_regul.genes[,1:3],up_regul.genes[,6],up_regul.genes[,5],up_regul.genes[,4])
	down_regul.genes<-anno.filter[down_regul.table$id,]
	down_regul.genes<-cbind(down_regul.genes[,1:3],down_regul.genes[,6],down_regul.genes[,5],down_regul.genes[,4])
	print(paste("Up regulated exon number is: ",dim(up_regul.table)[1]))
	print(paste("Down regulated exon number is: ",dim(down_regul.table)[1]))
	returnResult<-list(upRegulGenes=up_regul.genes, upRegulTables=up_regul.table, downRegulGenes=down_regul.genes, downRegulTables=down_regul.table)
	return(returnResult)
}
##NOISeq:
##e.g.resultByNOISeq<-sigExonByNOISeq(count.table.1stExon,anno=count.table.1stExon.GC.geneLen)
sigExonByNOISeq<-function(count.table,group=factor(c(1,1,2,2)),filterNum=1,fdr=0.05, anno=NULL){
	filter <- apply(count.table,1,function(x) sum(x)>=filterNum)
	common <- intersect(rownames(anno),rownames(count.table[filter,]))
	count.table.filter<-count.table[common,]
	anno.filter<-anno[common,]
	myfactors = data.frame(conditions = factor(c(1,1,2,2)))
	mychroms = data.frame(Chr = anno.filter[,1],GeneStart=anno.filter[,2],GeneEnd=anno.filter[,3])
	rownames(mychroms)=rownames(anno.filter)
	mylength=anno.filter[,3]-anno.filter[,2]
	names(mylength)=rownames(anno.filter)
	mydata <- readData(data = count.table.filter, length = mylength, biotype = NULL, chromosome = mychroms, factors = myfactors)
	mynoiseq = noiseq(mydata, k = 0.1, norm = "rpkm", factor = "conditions", replicates = "biological")
	res = mynoiseq@results[[1]]
	up_regul.table<-res[res$foldChange > 2 & res$padj < fdr,]
	down_regul.table<-res[res$foldChange < 0.05 & res$padj < fdr,]
	up_regul.genes<-anno.filter[rownames(up_regul.table),]
	up_regul.genes<-cbind(up_regul.genes[,1:3],up_regul.genes[,6],up_regul.genes[,5],up_regul.genes[,4])
	down_regul.genes<-anno.filter[rownames(down_regul.table),]
	down_regul.genes<-cbind(down_regul.genes[,1:3],down_regul.genes[,6],down_regul.genes[,5],down_regul.genes[,4])
	print(paste("Up regulated exon number is: ",dim(up_regul.table)[1]))
	print(paste("Down regulated exon number is: ",dim(down_regul.table)[1]))
	returnResult<-list(upRegulGenes=up_regul.genes, upRegulTables=up_regul.table, downRegulGenes=down_regul.genes, downRegulTables=down_regul.table)
	return(returnResult)
}


###make density scatter plot
library(fields)
fudgeit <- function(){
	xm <- get('xm', envir = parent.frame(1))
	ym <- get('ym', envir = parent.frame(1))
	z  <- get('dens', envir = parent.frame(1))
	colramp <- get('colramp', parent.frame(1))
	image.plot(xm,ym,z, col = colramp(256), legend.only = T, add =F)
}
Lab.palette <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

makeScatterDensityplot<-function(x,filename="test.pdf",xlim=c(0,100),ylim=c(0,100),nrpoints=0,...){
	pdf(filename, paper="special", height=4, width=4)
	par(plt = c(0.1171429 ,0.8200000, 0.1457143, 0.8828571))
	par(mar = c(5,4,4,5) + .1)
	plot(0,0,xlim=xlim,ylim=ylim,xlab="",ylab="")
	rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col ="#00007F")
	par(new=T,mar = c(5,4,4,5) + .1)
	smoothScatter(cbind(x[,1],x[,2]),nbin = 1000,bandwidth=0.001,colramp = Lab.palette,xlim=xlim,ylim=ylim, postPlotHook = fudgeit,nrpoints = nrpoints,...)
	abline(xlim[1],1)
	dev.off()
}

makeScatterDensityplotDefault<-function(x,filename="test.pdf",xlim=c(0,100),ylim=c(0,100),...){
	pdf(filename, paper="special", height=4, width=4)
	par(plt = c(0.1171429 ,0.8200000, 0.1457143, 0.8828571))
	par(mar = c(5,4,4,5) + .1)
	correlation<-paste("cor=",cor(x)[1,2],sep="")
	num<-paste("Num:",dim(x)[1],sep="")
	#smoothScatter(x,nbin = 1000,bandwidth=0.001,main=paste(correlation,num,sep="\n"),xlim=xlim,ylim=ylim, ...)
	smoothScatter(x,main=paste(correlation,num,sep="\n"),xlim=xlim,ylim=ylim, ...)
	abline(xlim[1],1)
	dev.off()
}

scatterPlotBetweenTwoBed<-function(x1,x2,outputNames="test.pdf",colToUse=4,...){
	x1.data<-read.table(x1,sep="\t",header=F,skip=1)
	x2.data<-read.table(x2,sep="\t",header=F,skip=1)
	rownames(x1.data)<-paste(x1.data[,1],x1.data[,2],x1.data[,3],sep=":")
	rownames(x2.data)<-paste(x2.data[,1],x2.data[,2],x2.data[,3],sep=":")
	common<-intersect(rownames(x1.data),rownames(x2.data))
	makeScatterDensityplotDefault(cbind(x1.data[commong,colToUse],x2.data[commong,colToUse]),filename=outputNames,...)
}

##make average plot for any two provided matrix(check consistency between replicates):
makeAveragePlot<-function(x1,x2,filename="test.pdf",scale=1000,step=20,bin_size_align=1,colors=c("black","purple"),line_types=c(1,1),legendNames=c("rep1","rep2")){
	dataSeq<-seq(((length(x1[1,]))/2)-as.integer(scale/bin_size_align), ((length(x1[1,]))/2)+as.integer(scale/bin_size_align), by=as.integer(step/bin_size_align))
	axisSeq<-seq(0-scale, scale, by=step)
	valueGch1<-NULL
	valueGch2<-NULL
	for(i in dataSeq){
		if(floor(i+step/(2*bin_size_align)-1)-ceiling(i-step/(2*bin_size_align)) <= 0){
			valueGch1<-cbind(valueGch1,mean(x1[,ceiling(i-step/(2*bin_size_align))], na.rm=T))	
			valueGch2<-cbind(valueGch2,mean(x2[,ceiling(i-step/(2*bin_size_align))], na.rm=T))	
		}else{
			valueGch1<-cbind(valueGch1,mean(colMeans(x1[,ceiling(i-step/(2*bin_size_align)):floor(i+step/(2*bin_size_align)-1)], na.rm=T), na.rm=T))
			valueGch2<-cbind(valueGch2,mean(colMeans(x2[,ceiling(i-step/(2*bin_size_align)):floor(i+step/(2*bin_size_align)-1)], na.rm=T), na.rm=T))
		}
		
	}
	pdf(filename, paper="special", height=4, width=4)
	par(oma=c(1, 1, 1, 1))
	par(mar=c(1, 1, 3, 1))
	plot(axisSeq,valueGch1,type="l",axes=FALSE,xlab="",ylab="",ylim=c(0,100),col=colors[1],lty=line_types[1],font=2,lwd=3)
	par(new=T)
	plot(axisSeq,valueGch2,type="l",axes=FALSE,xlab="",ylab="",ylim=c(0,100),col=colors[2],lty=line_types[2],font=2,lwd=3)
	axis(2,at=seq(0,100,by=20),labels=format(seq(0,100,by=20),digits=3),lty=1,font=2,cex.axis=1.0,cex.lab=1.0,font.lab=2,lwd=2)
	axis(1,at=seq(0-scale, scale, by=scale/2),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2)
	mainTitle=paste("numElemCenterToAlign:",length(x1[,1]))
	title(mainTitle, cex.main = 0.6, font.main= 4, col.main= "black",xlab="Distance to elements (bp)")
	legend("topright",legendNames, col=colors,lty=line_types,cex=0.5,lwd=2)
	
	abline(v=0)
	dev.off()
}


returnMatrixOnBed<-function(x1,x2){
	bed1<-read.table(x1,sep="\t",header=F,skip=1)
	bed2<-read.table(x2,sep="\t",header=F,skip=1)
	rownames(bed1)<-paste(bed1[,1],bed1[,2],bed1[,3],bed1[,6],sep=":")
	rownames(bed2)<-paste(bed2[,1],bed2[,2],bed2[,3],bed2[,6],sep=":")
	common<-intersect(rownames(bed1),rownames(bed2))
	value<-cbind(bed1[common,4],bed2[common,4],bed1[common,5],bed2[common,5])
	rownames(value)<-common
	value
}

returnMatrixOnBedGraph<-function(x1,x2){
	bed1<-read.table(x1,sep="\t",header=F,skip=1)
	bed2<-read.table(x2,sep="\t",header=F,skip=1)
	rownames(bed1)<-paste(bed1[,1],bed1[,2],bed1[,3],sep=":")
	rownames(bed2)<-paste(bed2[,1],bed2[,2],bed2[,3],sep=":")
	common<-intersect(rownames(bed1),rownames(bed2))
	value<-cbind(bed1[common,4],bed2[common,4])
	rownames(value)<-common
	value
}


#####choose the nearest gene's FPKM
selectNearestGenes<-function(x1,x2,x3,x4, numNearest=5){
	require(plyr)
	x1.value<-read.table(x1,sep="\t",header=F)
	rownames(x1.value)<-paste(x1.value[,1],x1.value[,2],x1.value[,3],x1.value[,7],x1.value[,8],x1.value[,9],x1.value[,11],sep=":")
	x2.value<-read.table(x2,sep="\t",header=F)
	rownames(x2.value)<-paste(x2.value[,1],x2.value[,2],x2.value[,3],x2.value[,7],x2.value[,8],x2.value[,9],x2.value[,11],sep=":")
	x3.value<-read.table(x3,sep="\t",header=F)
	rownames(x3.value)<-paste(x3.value[,1],x3.value[,2],x3.value[,3],x3.value[,7],x3.value[,8],x3.value[,9],x3.value[,11],sep=":")
	x4.value<-read.table(x4,sep="\t",header=F)
	rownames(x4.value)<-paste(x4.value[,1],x4.value[,2],x4.value[,3],x4.value[,7],x4.value[,8],x4.value[,9],x4.value[,11],sep=":")
	common<-intersect(rownames(x1.value),rownames(x2.value))
	common<-intersect(common,rownames(x3.value))
	common<-intersect(common,rownames(x4.value))
	x.merge<-cbind(x1.value[common,],x2.value[common,12],x3.value[common,12],x4.value[common,12])
	
	M<-cbind(abs(x.merge[,8]-x.merge[,2]),abs(x.merge[,8]-x.merge[,3]),abs(x.merge[,9]-x.merge[,2]),abs(x.merge[,9]-x.merge[,3]))
	x.merge<-cbind(x.merge,apply(M,1,min))
	x.merge<-x.merge[order(x.merge[,1],x.merge[,2],x.merge[,3],x.merge[,17]),]
	x.merge.data<-data.frame(chr=x.merge[,1],start=x.merge[,2],end=x.merge[,3],name=x.merge[,10],dist=x.merge[,17],HCT116_F=x.merge[,12],HCT116_A=x.merge[,14],DKO1_F=x.merge[,15],DKO1_A=x.merge[,16])
	s<-ddply(x.merge.data, .(chr, start, end), function(x) x[1:min(numNearest,nrow(x)),])
	s.HCT116_F.mean<-ddply(s, .(chr, start, end), function(x) mean(x[1:nrow(x),6]))
	s.HCT116_A.mean<-ddply(s, .(chr, start, end), function(x) mean(x[1:nrow(x),7]))
	s.DKO1_F.mean<-ddply(s, .(chr, start, end), function(x) mean(x[1:nrow(x),8]))
	s.DKO1_A.mean<-ddply(s, .(chr, start, end), function(x) mean(x[1:nrow(x),9]))
	s.mean<-cbind(s.HCT116_F.mean,s.HCT116_A.mean[,4],s.DKO1_F.mean[,4],s.DKO1_A.mean[,4])
	s.mean
}

makeBarChart3SamplesPlot<-function(x,name="test.barChart.pdf", colList=c("black","orange","red")){
	pdf(name, paper="special", height=5, width=5)
	overlapped<-x[,4:length(x[1,])]
	no_express_criteria=0.1
	low_express_criteria=c(0.1,1)
	high_express_criteria=c(1,max(overlapped))
	par(mar=c(4.5,5,6,0))
	names=rep("",4)
	totalData<-NULL
	for(i in c(1:4)){
		no_express<-100*length(overlapped[overlapped[,i]<=no_express_criteria,][,i])/length(overlapped[,i])
		low_express<-100*length(overlapped[overlapped[,i]>low_express_criteria[1] &overlapped[,i] <= low_express_criteria[2],][,i])/length(overlapped[,i])
		high_express<-100*length(overlapped[overlapped[,i]>high_express_criteria[1] &overlapped[,i] <= high_express_criteria[2],][,i])/length(overlapped[,i])
		totalData<-cbind(totalData,c(no_express,low_express,high_express))
	}
	barplot(totalData, names.arg=names,col=colList,xlab="Samples",ylab="Percentage of genes",main="",ylim=c(0,100),cex.axis=1.5,font.axis=2,cex.lab=2,font.lab=2,cex.main=2,font.main=2,horiz=F,space=c(0.5,0.1,0.5,0.1))
	dev.off()
}

###add error bar to barplot:
error.bar <- function(x, y, upper, lower=upper, length=0.1,...){
	if(length(x) != length(y) | length(y) !=length(lower) | length(lower) != length(upper))
		stop("vectors must be same length")
	arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
}
y <- rnorm(500, mean=1)
y <- matrix(y,100,5)
y.means <- apply(y,2,mean)
y.sd <- apply(y,2,sd)
barx <- barplot(y.means, names.arg=1:5,ylim=c(0,1.5), col="blue", axis.lty=1, xlab="Replicates", ylab="Value (arbitrary units)")
error.bar(barx,y.means, 1.96*y.sd/10)



#######heatmap utils



